Lab 6 - Network Optimization Models - The Minimum Spanning Tree Problem

Information on group members:

1. 156035, Kuba Czech
   
2. 156045, Wojciech Nagórka

%matplotlib inline
import csv
import matplotlib.pyplot as plt
import numpy as np

1.1) Get acquainted with the below "City" class definition. It has three major fields: name, city's longitude, and latitude. The coordinate's input is in the format "00°00'E/N," and it is converted to a list [degrees, minutes]. Also, these coordinates are also simply converted and stored as x and y coordinates.

class City:
    def __init__(self, name, long, lat):
        self.name = name
        self.long = [float(long.split("°")[0]), float(long.split("°")[1].replace("'",'').replace("E",''))]
        self.x = self.long[0]*60 + self.long[1]
        self.lat = [float(lat.split("°")[0]), float(lat.split("°")[1].replace("'",'').replace("N",''))]  
        self.y = self.lat[0]*60 + self.lat[1]
       
    def __repr__(self):
        return str(self.name) + " " + str(self.lat) + " " + str(self.long)

1.2) The below piece of code loads cities data from the cities.csv file and creates a list of cities.

cities = []

with open('cities.csv', newline='', encoding="utf-8") as csvfile:
    r = csv.reader(csvfile, delimiter=',', quotechar='|')
    for row in r:
        cities.append(City(row[0], row[1], row[2]))
        
for c in cities[:10]: ### TEST
    print(c)

Adamów (siedleckie) [51.0, 45.0] [22.0, 15.0]
Adamów (zamojskie) [50.0, 36.0] [23.0, 10.0]
Adamówka [50.0, 16.0] [22.0, 42.0]
Aleksandrów [51.0, 16.0] [19.0, 59.0]
Aleksandrów Kujawski [52.0, 53.0] [18.0, 42.0]
Aleksandrów Łódzki [51.0, 49.0] [19.0, 19.0]
Alwernia [50.0, 4.0] [19.0, 32.0]
Andrespol [51.0, 44.0] [19.0, 37.0]
Andrychów [49.0, 52.0] [19.0, 20.0]
Andrzejewo [52.0, 50.0] [22.0, 12.0]

1.3) The function below plots all cities according to their X and Y coordinates. Note that there are so many cities that Poland's boundaries can be observed. Feel free to modify the code. The unrequired input "edges" is a list of pairs of indices indicating edges to be drawn.

def plotMap(cities, edges=[]):
    fig=plt.figure(figsize=(6,6), dpi= 100, facecolor='w', edgecolor='k')
    
    for i, e in enumerate(edges):
        color="red"
        X = [cities[e[0]].x, cities[e[1]].x]
        Y = [cities[e[0]].y, cities[e[1]].y]
        plt.plot(X, Y, lw=1, ls="-", marker="", color=color)
        
    X = [c.x for c in cities]
    Y = [c.y for c in cities]
    plt.plot(X, Y, ms=2, ls="", marker="o")
    
    plt.show()

plotMap(cities, edges=[[1,200]])

1.3) The dataset consists of 2249 cities. It may be too much when it comes to implementing the algorithm for constructing the minimum spanning tree (computational burden). Therefore, you can consider only every n-th row from the dataset to ensure that the computations will be completed in a reasonable time.

n = 3
cities_small = cities[2::n]
plotMap(cities_small, edges=[[0,200]])

1.4) Finish the function below for constructing a distance matrix. Each i-th row should contain distances from the i-th city to other cities. For convenience and to improve the algorithm's computational complexity for constructing the spanning tree, it is suggested to keep rows sorted according to distance. Rows may contain information on the j-th indices of the respective destinations along with the distances in a tuple (j-th indice, distance).

You can use the Euclidean distance for computations. You can use the x and .y fields that store the geographical coordinates. This way of formulating a distance does not have much sense because the computation would be based on angles. Nonetheless, it is allowed in this exercise to do so, but if you wish, you can calculate the distances correctly by computing arc lengths between two points on a sphere - Earth.

def getSortedDistances(cities):
    D = []

    for i, city1 in enumerate(cities):
        l = []
        for i, city2 in enumerate(cities):
            # longituteds and latitudes are in degreees
            lat1 = np.deg2rad(city1.lat[0] + city1.lat[1]/60)
            lat2 = np.deg2rad(city2.lat[0] + city2.lat[1]/60)

            long1 = np.deg2rad(city1.long[0] + city1.long[1]/60)
            long2 = np.deg2rad(city2.long[0] + city2.long[1]/60)

            dist = np.arccos((np.sin(lat1) * np.sin(lat2)) + (np.cos(lat1) * np.cos(lat2) * np.cos(long1 - long2))) * 6371

            l.append((i, dist))
        l = [i for i in l if i[1] != 0.0]
        l.sort(key=lambda x: x[1])
        D.append(l)
    return D

D = getSortedDistances(cities_small)
print(D[0][:5])

/var/folders/8m/nx_b_wh17dg77b9kb95gxng80000gp/T/ipykernel_44983/62982253.py:14: RuntimeWarning: invalid value encountered in arccos
  dist = np.arccos((np.sin(lat1) * np.sin(lat2)) + (np.cos(lat1) * np.cos(lat2) * np.cos(long1 - long2))) * 6371

[(550, 11.00036580940852), (287, 16.84614692625338), (27, 17.337232288265085), (474, 17.55795070311418), (300, 18.95399028168302)]

1.5) Auxiliary function for finding a corresponding index in the data set for a provided city name.

def getIdx(name, cities):
    for i, c in enumerate(cities):
        if name == c.name: return i
    return None

print(getIdx("Hel", cities_small))

178

1.6) Complete the below greedy algorithm for constructing the minimum spanning tree. As an input, it should return a list of edges of the constructed subgraph, which can then be used in the plotMap function to show the tree structure. You can set the city of Hel as an initial node.

def getSpanningTree(stratingNode, cities, D):
    N = len(cities)
    IN_TREE = set([stratingNode])
    edges = [[] for i in range(N - 1)]

    for _ in range(N-1):
        nodes_to_choose = []
        for i in IN_TREE:
            for idx, dist in D[i]: # i is in tree and idx is node to be added
            # 1. In each row find a node closest to the tree
                if idx not in IN_TREE:
                    nodes_to_choose.append([i, idx, dist])
                    break
        # 2. From candidates choose the best one
        nodes_to_choose.sort(key = lambda x: x[2])
        i, idx, dist = nodes_to_choose[0]
        edges[len(IN_TREE) - 1].extend([i, idx])
        IN_TREE.add(idx)

    return edges

edges = getSpanningTree(getIdx("Hel", cities_small), cities_small, D)
plotMap(cities_small, edges=edges)

1.6) Modify the code so that it displays the results after 10 intermediate steps. In this way, you can get a better insight into how the algorithm constructs the spanning tree. Additionally, you can use, e.g., a gradient to distinguish between different iterations of the algorithm when plotting the map (e.g., 1-st iteration = 100% red, last iteration == 100% green).

def getSpanningTreeFirstM(stratingNode, cities, D):
    N = len(cities)
    IN_TREE = set([stratingNode])
    edges = [[] for i in range(len(cities) - 1)]

    for _ in range(len(cities) - 1):
        nodes_to_choose = []
        for i in IN_TREE:
            for idx, dist in D[i]:
                if idx not in IN_TREE:
                    nodes_to_choose.append([i, idx, dist])
                    break
        nodes_to_choose.sort(key = lambda x: x[2])
        i, idx, dist = nodes_to_choose[0]
        edges[len(IN_TREE) - 1].extend([i, idx])
        IN_TREE.add(idx)

        if (len(IN_TREE) - 1) % 10 == 0 or len(IN_TREE) == len(cities):
            plotMap(cities, edges = edges[:len(IN_TREE) - 1])

getSpanningTreeFirstM(getIdx("Hel", cities_small), cities_small, D)